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ABSTRACT 





In this study, several turbulence models are analyzed 








in a 3-D finite element model of a micro-heat exchanger. 
The micro-heat exchanger consists of a narrow planar flow 


passage between parallel walls with small cylindrical pin 








fins spanning these walls with axes perpendicular to the 


direction of §$ flow. Turbulence model performance is 





compared with baselin xperimental data available in the 
literature that cover a range of low turbulent Reynolds 
numbers and spacing configurations. The metric for these 
comparisons is an array averaged Nusselt Number. 
Adjustments made to the coefficients in the turbulence 
models are explained in terms of their physical 
significance to the complex flow environment of a pin fin, 
cross flow, micro-heat exchanger. Applications of this 
research include cooling of turbine blades and of closely 


spaced electronics. 
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L3 INTRODUCTION 


High temperatures on power turbine blades have been a 
limiting factor in advancing the design of gas turbine 


engines. Conventional cooling schemes in advanced gas 





turbine engines will be unable to remov th heat 








transferred from the mass stream fast enough to prevent 





damage to the blades. A new cooling scheme has been 


proposed which may be adaptable to the turbine blade 





cooling problem. It consists of a narrow flow passage 


between two endwalls is separated by rows of short pin fins 








with axes perpendicular to the coolant flow. It is 


envisioned that this flow passage would be tightly 





"wrapped" around a turbine blade forming a_ protective 
cooling shroud. Bleed-off air for cooling would be 


supplied through the blade root and distributed at the 








leading edge of the passage flow with an exit to the 


turbine flow at the trailing edge of the blade. This paper 








addresses the turbulence models used in the 3-D numerical 
modeling of two arrangements of this cooling scheme over 
several low, but fully turbulent Reynolds numbers. The 


turbulence models have been fine tuned to agree with 








availabl xperimental data in the literature, and optimal 





models and coefficients have been identified. 
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II. BACKGROUND AND OBJECTIVES 


A. HEAT EXCHANGER CONSTRUCTION 


The micro-heat exchanger modeled in this study is made 
up of a narrow flow passage supported by 10 staggered rows 
of pin fins perpendicular to the flow. Figure 1 is a 
definition sketch of a typical configuration, similar to 


the type constructed by Metzger, et al. (1982) which forms 





the basis of corroboration for the current numerical study. 
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0 


W 


— 


FLOW 


> 





D 





Figure 1. Definition sketch of a staggered pin- 
fin array 


The pins are arranged with a constant axial pitch (xX, 


distance between pin centers in the flow direction) to form 











the reference length scale. The axial pitch was chosen to 
be 12.7mm to be consistent with earlier studies. Two 


different diameter pins were used in this setup, 8.46mm and 


3 








5.08mm. The X/D and S/D ratios defin th relative 





streamwise and spanwise pin spacing respectively for a 
general model. As a result of these settings the 
transverse pitch (S) was not a constant from model to 


model, although axial pitch was. 


As in the work of Metzger et al (1982), the endwalls 





(the broad upper and lower plates) were modeled as 


isothermal surfaces. The pins were modeled as isothermal 





at the temperature of the bounding endwalls as a first 














approximation. It can be shown with a standard pin fin 
relation (Incropera and Dewitt, 1996) that the isothermal 
assumption for the pins is a valid assumption. The side 
walls were modeled as insulated, again to duplicate 





Metzger's work and an insulated inlet and outlet section 





were added (not shown in Figure 1). 


B. PREVIOUS WORK 


The staggered pin-fin array geometry, in one form or 


other, serves as the basis for a large number of compact 





heat exchanger configurations, and has received 








considerable attention in the past. Th xperimental work 


of VanFossen (1982), Metzger et al (1982), Chyu (1990), 





Chyu and Goldstein (1991), and Jubran et al. (1993) and 





others have explored various issues such as the differences 


between short and long pins, the contribution of pin heat 





transfer versus endwall heat transfer, the presence of an 


optimal streamwise and spanwise spacing/pitch, etc. 


Only recently, with the advent of more powerful 





computational machines and software tools, have numerical 


studies of compact heat exchangers been attempted as 


reviewed by Shah et al. (2001). For the staggered pin fin 








array configuration in particular, Donahoo et al. (2001) 


have used a general purpose viscous flow solver to simulate 





the fluid flow and heat transfer behavior in a 2-D 











numerical study. However the 2-D nature cannot be used to 


treat the endwalls where the pins meet the surface. 


Th xperimental baseline for this study is the 





Metzger et al. work of 1982. Metzger showed the variation 
of streamwise heat transfer, overall array heat transfer, 
and overall flow friction in staggered, short pin fin 
arrays. His experimental work has provided the baseline 


for this computer model. His experimental work showed that 








there was a substantial difference between short pin fin 





arrays and long cylinder (tube bank) arrays. Long 





cylinders had been the previous model for the array. The 








relations Metzger developed for overall heat transfer 


showed a lower array averaged Nusselt number over the array 





based on his experiments. The long cylinder predictions 





did not correctly predict the convection coefficients 


observed in the short pin fin array. 


Cc. OBJECTIVE 


The primary objective of this study was to select a 





turbulence model that accurately reflects th xperimental 





results seen in previous work. As part of the selection 


process, the turbulence model had to be optimized since the 





default coefficients in a model are typically better suited 


for flat plate or long cylinder tube-bank type open 





geometries, but not very applicable to confined tortuous 


flow geometries such as the short pin fin array in this 


5 


study. Another objective is to leave the code open for 
future work where pin shape and spacing/layout could be 
changed. The goal would be to arrive at a recommendation 
for the most optimal turbulence model/s and their 


coefficients that would provide the best corroboration with 








availabl xperimental data. Furthermore it is hoped that 








subject to futur xperimental work, the current findings 
can be extended to have validity over a wide range of 


Reynolds numbers and geometrical configurations thus 





providing greater applicability and confidence in future 





numerical simulations. 


ue i TURBULENCE MODELS 


A. INITIAL MODELING 


1. Boundary Conditions 





The numerical model of Metzger's xperiment consists 


of two different boundaries. Fluid elements wer defined 





which had rigid boundaries, eliminating the need _ for 
defining a rigid structure to surround the fluid elements. 


To reduce the number of nodes required for modeling, 





symmetry was exploited. A horizontal plane, parallel to 
the flow direction, halfway between the upper and lower 
isothermal boundaries of the passage defines an adiabatic 
boundary. No net properties cross this boundary. The 
length scale in the flow passage is sufficiently small that 
buoyancy effects are negligible. This boundary was defined 


as insulated and velocity in the vertical direction was 





Zero. However the no-slip condition along the plane was 


released. A similar symmetry plane was defined 





perpendicular to the first symmetry plane and parallel to 





the flow direction. The boundary was again modeled as 
insulated, crossing velocity was set to zero and the no- 


slip condition along the plane was released. These two 





planes allowed a model to only have nodes in one quadrant 





as the other quadrants would produce the same result. 





Observing the downstream direction, the lower-right 


quadrant was the guadrant chosen for modeling. 





Symmetry 
Planes 








Inle 


Endwall Us Passage Wall 





Figure 2. Modeled Portion of a Pin Fin Heat 
Exchanger 





12.7 cm xI(2D) Pin-Fin Array (10 rows) xa2D) 12.5 cm 


Figure 3. Actual model used for experiment 





The pins were modeled as isothermal and so was the 


supporting endwall. The remaining wall, the sidewall, 





parallel to the pin axes was insulated. The no slip 
condition was applied to all boundaries except the symmetry 
planes. Straight duct entrance and exit regions were 


attached at the inlet and outlet sections of the array test 





section to be consistent with the experimental rig of 





Metzger et al (1982). They were modeled slightly longer 


than the Metzger model to prevent the outlet conditions 





from affecting the test section. 


The isothermal boundaries were fixed at 306 K. Dry 





air (Pr = 0.7) was used as working fluid, entering the 


k, 
inlet section at 300K, p =1.1769 28. , u=1.8468x10°—2 , 
m m:s 


K =2.6294x10?——, and C, =10040——. Inlet velocity was 


m:K 2 





specified at the entrance to the test section based on a 


maximum velocity type Reynolds number. 


The maximum velocity Reynolds number is defined as 


pec Oe! (Metzger et al., 1982), where V,, is the 
U 

inviscid, geometry based maximum velocity the fluid reaches 

when passing through the test section. D is the pin 

diameter. 


2. Initial Conditions 


Bach run was initiated with the walls and pins of the 


test section at 306K and the incoming air at 300K. Inlet 








velocity was determined by th specified Reynolds number 








and fluid properties. 


3. Mesh Structure 


The mesh was specified to concentrate nodes in the no 
slip boundary areas of the model, where flow and 
temperature gradients might be expected to be highest (Fig 
4). Around the pins, nodes were concentrated in the high 


velocity regions where conditions change the most (Fig 5). 





71 6 2002 
12:03:53 


Re=7310, $/D=2.5,X/D=1,5, To =306, h/v/pH5/16/2,COLG, 250,3KE 





Figure 4. Standard Mesh showing nodes 
concentrated in the areas of the no slip 
boundary. 


Re=7310, $/D=2.5,X/D=1,5, To =306, h/v/pHa5/ 17/2, CLG, 250,3KE 





Figure 5. Standard Mesh (Top View) with nodes 
concentrated around the high velocity regions of 
the pins 
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Nodes were also concentrated along the endwall. 


For constructing the numerical model of the heat 
exchanger and for the majority of the data gathering, the 


finite element based modeling and simulation package, ANSYS 





was used. The versions available during this research were 
ANSYS version 5.7, 6.0, and 6.1. The bulk of the research 
was done using ANSYS 6.0. The ANSYS license at the Naval 
Postgraduate School (NPS) limited models to 256,000 nodes. 








Several runs were attempted using the Department of Defense 


High Performance Computing Major Shared Resource Center 





(HPCMSRC). The ANSYS license at the HPCMSRC was node 


unlimited. 








At the entrance and exit of the test section, as well 


as midway between each row of pins (apparent in figures 4 





and 5), a vertical plane of nodes was constructed. This 





allowed bulk temperature and mass flows to be calculated 
independent of the ANSYS program. Such a layer of nodes 
would also facilitate future work where pin shape and 


layout are expected to be changed. 


B. SKE MODEL 


The Standard k-€ (SKE) Model is the default turbulence 


model for ANSYS 6.0. The SKE model is a Reynolds Averaged 





Navier-Stokes (RANS) type of Eddy Viscosity Model (EVM). The 
SKE model uses 2 equations to close the governing 
equations. The EVM approach introduces a turbulent 


viscosity term, U,, which relates the fluctuating velocities 


of a turbulent flow to a viscous stress term, 


ah 


. Uu. 
—uu,=——, in tensor notation, which is used in the 


Reynolds stress terms of the time averaged, turbulent, 





incompressible Navier-Stokes gquations, as well as in the 
energy equation and the continuity equation. Asan 


example, the x component of the Reynolds stress term is: 


Ofoaren\ tae. 0. one ee: 
of =-<( vv. )-—(pv.v,]-—(ev.v.], where V.,V.,V. are the 

ox oy dz 
fluctuating components of the velocities. The Reynolds 


stress terms contain the indeterminate part of the Navier- 





Stokes relations, the mean of the product of the 
fluctuating velocity components. Replacing that product 


with the relation involving the turbulent viscosity allows 








the Navier-Stokes momentum equations, continuity and energy 





equations to be reshaped. In the SKE model two transport 


equations are derived from the governing equations, the 





Turbulent Kinetic Energy equation: 


apk A(pV.k)  A(pV,k) | A(pv.k) 
ot ox oy oz 


3b bok acd Lok foes Mok +u,® 
dx|o, dx} dylo, day } dazlo, dz Aaa 


“pe ors =| 


+g,=+g. 
* Ox &y oy aa 


t 





and the Dissipation Rate equation: 


ape , 2(P%e)  A(pV.e) | A(pve) 


ot ox oy a 
d(u,de) dfu,de) Af u, dE € e? 
Be Pt Schco omc Wnt cai ipa aed 2c) oy MT We eee ga, Cam 
se sae | op (Soe | ae eee | He 2P k 
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+ 


Bia §y 3 825 
o ox” oy dz 


t 


J- 
C, ( ama or, or, =) ne 


where k is the turbulent kinetic energy, Mra Noa Vb ee 


Re ‘Instantaneous (non-fluctuating). velocities, <6, is: the 





Schmidt number for the transport of turbulent kinetic 


energy, rc is the viscous dissipation term 


du. ou, \du, ; 
nu{ pees —|, € is the turbulent kinetic energy 





dissipation rate, C, and C, are the buoyancy constants, 68 





is the coefficient of thermal expansion, oO, is the turbulent 


t 


Prandtl number, g is the acceleration scalar, JT is the 





static temperature, oO is the Schmidt number for the 


€ 





transport of turbulent kinetic energy dissipation rate, C 





le 


is the shear rate generation term multiplier, C, is scaled 


dittusivity rate multiplier, C, is. the: turbulent. viscosity 


constant, and WU, is the turbulent viscosity broken down in 


2 
Table 1 contains the values of 





k 
the equation to HW,=C P 
E 


iv 





the SKE constants: 


1.3 




















Table 1. Constants for SKE Model 


Solving these equations at each node in the model for 








a set number of iterations yields values for k and € which 





are used to calculate wu, from the previous equation. U, is 





used to calculate the indeterminate fluctuating velocities 





in the Navier-Stokes quations. The Navier-Stokes 


equations are then solved for a new set of mean velocities. 


This new set of velocities, with fluctuating 
components, is then used to solve for temperature in a 


first law relation: 














c| ——+V_—4+V —+V. = at 
f : ax? dy? dz? 


OU OF ap Oh cg OE to Of 01 OL ei 
ot a& dy *a 


14 


where K is the thermal conductivity of the fluid, 
assumed to be isotropic, c is the specific heat, and the 


velocities ViVeN, contain the mean and fluctuating 


components. 


The local heat flux, ds can be estimated at non- 


insulated or adiabatic boundaries by the following relation 


where y is the direction normal to the transfer surface, 





oT 
q =-K—  . From this, Newton's Law of Cooling leads to a 
oy y=0 
convection coefficient, a where T, is the surface 
(7, -T, ) 
temperature and T, is the fr stream temperature. Since 
the flow complex, ANSYS 6.0 uses T, in place of Tf. fA 
the bulk temperature of the adjacent fluid element. 





The ultimate goal of the simulation is to calculate an 





array averaged Nusselt number, so an overall array 





convection coefficient, h, is calculated based on the total 


heat transferred to the fluid. 


q 


=—— (1.4) 
A:LMTD 


where q is the total amount of heat in watts 


transferred across the boundaries, A is the wetted area of 








the test section, and LMTD is the Log-Mean Temperature 





Difference. 


LMTD =—————_ _ (1.5) 
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where TT, is the bulk inlet temperature, AT is 


calculated from APS. T is the constant isothermal 


mn 7 C wall 
P 





wall temperature for both endwalls and pins, and 


T,=1,+AT. The array averaged Nu is calculated from: 
Nu=— (1.6) 
K 


where D is the diameter of a pin and K is the thermal 
Ve ; io ere 
conductivity of air evaluated at <a a This formulation 


of Nusselt number is the most appropriate for comparing the 





experimental results of Metzger with the current ANSYS 


simulations. 


One drawback to the SKE model, as well as all the 2 


equation EVM's, is that it does not accurately predict 











turbulent kinetic energy and dissipation rate at elements 








along the wall. In this case a wall turbulence model is 





used. The ANSYS default model for wall turbulence is the 








Van Driest model. The Van Driest model is used at the wall 











nodes to determine viscosity. It assumes an eddy viscosity 





of zero at the wall and gives an accurate value for wall 











conductivity. ANSYS 6.0 uses the continuous Van Driest 








formulation with the velocity at a certain distance from 





the wall as a reference velocity. This form of the log law 


of the wall is solved for a wall viscosity from 
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where T is the shear stress, k is the slope of the 





wall parameter, FE is the Law of the Wall constant, U, is 


the laminar absolute viscosity, and 6 is the distance from 








the wall. The wall lement ffectiv VISCOSBTEY;. [oy ES 


derived from: 


Cc. RENORMALIZED GROUP (RNG) MODEL —- DIFFERENCES FROM SKE 


The RNG model uses an advanced statistical technique 


called Renormalization Group Theory applied to the 





instantaneous Navier-Stokes quations. The two transport 


equations derived for the SKE model appear in the RNG model 





with the exception of the C, multiplier for the kinetic 





energy generation term in the dissipation rate equation. 





In the RNG model the C, term takes the form 


C_ =1.42- = (1.9) 
le 1+ Bn 


where 1,=4.38, B=0.012, and n is a function of the 


stress deformation tensor, Si The tensor, Sir is a 


function of velocities in the flowfield, 
Lf av, eV, 


= —|——t+ Mes lO)) 
Dove OV 





There are differences in the buoyancy terms of the 





equations but for this heat exchanger as noted earlier 





buoyancy effects can be ignored. In the RNG model several 
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constants have different numerical values, the default 





numerical constants are summarized in Table 2: 





Table 2. RNG Model Coefficients 


This gives the RNG model an added level of complexity 
and capability over the SKE model due to the more 


distributed nature of accounting for the turbulent kinetic 








energy generation term in the flow field. In the SKE 








model, the kinetic energy generation coefficient in the 





dissipation rate equation is the same throughout the 





flowfield. In the RNG model, the coefficient is dependent 


upon the velocity components at each point in the flow 





field. ANSYS 6.0 uses a slightly more simplified version 





than that presented in Numerical Analysis of Compact Heat 


Exchanger Surfaces (R. K. Shah et al., 2001). 
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Past research, as noted R. K. Shah et al., describes 
RNG as being more accurate and reliable than other 2 


equation EVM models. Because it is a new model, relatively 








few studies have been done with it, but it appears to do 


better than SKE with low Reynolds numbers but worse when 





vortex shedding becomes an issue, such as with bluff bodies 
or flow around cylinders as in external flows with open 


geometries. Although the comparisons were based on 2D 





models only, RNG'sS success appears to be highly geometry 
dependent. 
D. NKE —- DIFFERENCES FROM SKE 





The New k-€ Model (NKE), also known as the Realizable 
k-€ Model, was proposed by T. H. Shih et al. in 1994. The 
model was primarily designed for high Reynolds number flows 
but has had some success for lower Reynolds number flows as 


well. The significant differences for a micro heat 





exchanger type of flow environment are that the NKE model 


uses a different formulation for the C, constant and a 


different form of the kinetic energy dissipation rate 





equation. 





The turbulent viscosity UW, is still calculated the sam 





way but the C, term is now a function of deformation 


tensors and antisymmetric rotation tensors. ANSYS 6.0 


replaces C, with this function: 





2 2 


C HAF LS ala, HO Ores | 4) pe (1.11) 


where the CQ ni term involves the angular velocity of the 





coordinate system which for this application is not 
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rotating, so that term drops out. Th ntire term under 





the radical is the mean velocity in the law of the wall 


model. 


The dissipation rate equation now becomes: 





ape  A(pvie) (pve) 9(pV.e) 
ot ox oy a 


0 / u, de d{ Uw, 0€ 0 { U, 0€ 1 
ete | hee PERG Nan kt ee ceeds uae ba, C oe 
aG ox } oy [ dy } 2(# Oz }re lee x2. 


2 ¢C, (1-C k 
a?  C,( Jer oro =) aes 











=C.p——+ | 
2P § Ax gy oy &. A 


k oO, 


and the C, multiplier now becomes: 


c,, =m c a. | (1.13) 











Ment 
where, 
k fav, ov, \ av. ov. 
=— || — + | + 1.14 
18 Fy vats aa aa 


and ¢.,. 1s 2 constant. 


The numerical values for the default NKE model are 


summarized in Table 3: 




















Table 3. NKE Model Coefficients 


E. FAILURES - ZEM, GIR, SZL 


There were three other turbulence models that were 


evaluated in this study, a zero equation model (ZEM), the 





non-linear model of Girimaji (GIR), and the Shih, Zhu, 











Lumley model (SZL). In the Zero Equation Model, the 
definition of a turbulent viscosity does not introduce 


another differential equation, and ANSYS calculates a 





characteristic length applied to the entire model. For 





this model ANSYS applied a default characteristic length of 
1 cm. Neither this characteristic length nor any changes 
made improved upon the performance of the default SKE model 


so it was removed from further consideration. 


The GIR model was more stable than the ZEM model but 





it provided results significantly different than the 








experimental data. It relies on a complex formulation for 


the C, term. 








The SZL model is a simpler model than NKE or GIR. t 








estimates a low level of turbulence for the flowfield and 





low effective velocity that can quickly become unstable in 


a complex flowfield. The C, term is significantly simpler. 











It was found that the heat transfer did not steady out 
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smoothly and it consistently overestimated the parameters 


used to determine the convection coefficient, h. 





Since each of these 3 models was found wanting 
compared to the earlier 3, they were not pursued in much 
further detail in this study. 

F. ADVECTION SOLVERS 


For this simulation the ANSYS recommendations for 





solving complex flow problems were used. Single solvers by 





themselves did not correctly simulate the flows’ but 


combinations of solvers did. For the first 30 iterations 





the Monotone Streamline Upwind (MSU) advection solver was 





used for momentum, turbulence, pressure, and temperature. 


MSU produces first order accurate solutions and is 





generally very robust. In the current study in particular, 








after the first 30 iterations, density was permitted to 
vary and the MSU was still used for momentum, turbulence, 
pressure, and temperature. For iterations 60 through 90 
the momentum was judged stable enough to switch to the 


Collocated Galerkin (COLG) method, an exact advection 





solver. For the iterations beyond 90, the COLG solver was 








used in the solution of momentum, turbulence, pressure, and 





temperature quantities for best consistency and stability 


in the final results. 
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IV. RESULTS AND DISCUSSION 


A. NUMERICAL EXPERIMENTS CONDUCTED 





In order for a typical numerical run to be considered 
successful, which could then be compared to other runs with 
different parameters, a number of numerical performance 
criteria had to be met. The first and most important was 
grid independence. To achieve this, the number of elements 
in the model was increased (in steps of at least 24000 


elements) until a less than two percent change was observed 





in the total heat transfer rate between two consecutive 


runs. Only one Reynolds number resulted in the heat 





transferred being outside the 2% limit, and that too was 


less than 3%. 


A second criterion was that no nodal temperature ever 


went above 306.4K in the last iteration in order to be 








consistent with the imposed boundary conditions. A third 
was that the ANSYS generated output bulk temp was within 
O.1K of the calculated temperature based on the mass flow 


rate. The final criterion that had to be met was that the 











sum of the differences in nodal temperatures between two 


consecutive iterations divided by the sum of the final 





nodal temperatures be less than 1x10" to ensure 


satisfactory convergence. 


Only when all these criteria were met was a run deemed 





complete, and hence ready for comparison with other 








completed runs. 


B. TURBULENCE MODEL PERFORMANCE 


Three models (SKE, RNG, NKE) that were able to meet 


the above robust benchmark criteria were compared. As 
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noted earlier, three others (ZEM, GIR, SZL) were also 
considered initially but failed to meet the basic benchmark 
criteria and were hence not pursued further. 


1. The SKE Model 


SKE is the default turbulence model implemented in 





ANSYS. The model generated stable solutions and in order 
to compare it with Metzger's experimental work an array 
averaged Nusselt number was computed for each variation of 
each model, and compared with its corresponding Metzger 


setup. Baseline results using SKE as the turbulence model 





consistently overestimat the Metzger Nusselt number with 








the exception of the high Reynolds number run in the 








tightest (X/D=1.5) configuration considered. Various 


coefficients in the SKE turbulent kinetic energy and 





dissipation rate eguations (1.1-.2)could be varied. Four 





out of a possible nine coefficients were chosen as the ones 


effecting the most critical changes while the others either 








dealt with buoyancy terms or were already accounted for by 








varying the chosen four. Each parameter was varied by an 


arbitrarily chosen value of 30% up and/or down from its 





default value and a run was made with only that parameter 


changed, i.e. all else was held fixed. After every run 





benchmarks wer checked and a comparison with Metzger's 


Nusselt number was made. Table 4 summarizes SKE 


performance with respect to changes in C, « 
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[nuns | x0 | so | anges | anges) mu | pre 
Number | X/D | S/D | Changed ao = aera 
| 3384 | 3228 | 483 | 
papas ef | a 
3980 15 25 C1-epsilon 39.83 32.28 23,39 
| 7310 | 1.5 | 25 | Default | 0 | 5058 | 4961 | 1.96 | 
2.5 |C1-epsilon -30 57.66 49.61 16.23 


P3600 [15 [25 Detaut [0 [| mim | 7775 | 77 


13800 | 1.5 | 2.5 |C1-epsilon| —-30 78.35 77.75 0.77 
| 3590 | 25 | 25 | Default | 0 | 2992 | 26.73 | 11.93 | 
3590 | 25 | 25 |Ct-epsilon| —_-30 36.13 26.73 | 35.17 
| 7340 | 25 | 25 | Default | 0 | 4569 | 4499 | 156 | 
7340 | 2.5 | 2.5 |Ct-epsilon| _-30 52.65 44.99 17.03 
| 13900 | 25 | 25 | Default | 0 | 70060 | 7161 | 1.41 | 


13900 2.5 2.5 |C1-epsilon +30 61.65 71.61 -13.91 
13900 2.5 2.5 |C1-epsilon -30 76.49 71.61 6.81 





Table 4. Summary of SKE performance when C, is 


varied 


The C, constant precedes the viscous dissipation term 





in the dissipation rate equation (1.2). Increasing the 











coefficient generally results in a lower approximation of 


the Nusselt number than the default case. Specifically the 





increase in C contributes to an increase in the source 


le 


term of the dissipation rate transport equation and results 








in a decrease in the turbulent kinetic energy. The case of 





Reynolds number of 13800 does not follow the same trend as 








the other Reynolds numbers. Figure 6 shows this decrease 


in turbulent kinetic energy as displayed by ANSYS 6.0, the 





view is from the adiabatic plane looking down towards the 





lower isothermal plane. The top figure is the default case 


25 


and the lower figure is with the 


increased C,, 


(note the scale of the colors). 


Figure 6. 


MOOAL SOLUTION 
STE?=4 
sim «1 

Java) 
RSf390 
Of «.7288-04 
OK «9,955 


» TeSE-D9 


1,106 3, 3916 


es 3960 ,5/052.5,4/061.5, eo =106,h/v/ p=5/1 7/2, COLG, SRE 


WOOAL SOLUTION 


STR}«4 
SUB #1 
pa +34 
RSS=0 
Sm =~. 1128-03 
SMC 6.04) 


[a7 


-i12£-03 


671319 e014 


PH 3950 5/052. 5, X/ PKL. 5, ho =306, hy pH Se 17/2, COG, SRE, C lepst30 


404 


1.343 2.685 


5.69 





402) 
4.356 








multiplier 


Changes in the turbulent kinetic energy 
in the simulated flow field due to a C,,. increase 


of 30% 


at Re=3980 
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In exploring: the. role -of “the -C, “coerfictent, At must 





be noted that C. also contributes to the source term of the 





dissipation rate transport equation (1.2) in a nonlinear 


manner and it was found that a decrease in oF leads to a 


decrease in turbulent kinetic energy at each node. The 
trend seen in; Table .5 shows: “that, decreasing .C,. leads -to 
lower Nusselt numbers than in the default SKE case, again 
with the exception of the case of Re = 13800. The effect of 


C, on the model increases at higher Reynolds numbers. 





Decreasing the turbulent Schmidt number for the 








transport of dissipation rate, O., affects only the 
diffusive terms of the dissipation rate transport equation 
(2) s A Gecredse: 18. @, leads to: ah: Increase» im the 
effective diffusivity of the dissipation rate. Table 5 
shows..that-.a- change “an oO, has wery .Jittlée impact -on- the 
default SKE Nusselt number. The model at Re=13900 was run 


on a different mesh than the default SKE model therefore no 





conclusions can be drawn from the large difference. 


The turbulent Schmidt number for the transport of 


kinetic energy, oO is the only parameter that affects the 


kl 


turbulent kinetic energy transport equation (1.1) through 





the diffusive terms. A decrease in Oo increases the 


k 
effective diffusivity of the turbulent kinetic energy. 


Table 5 shows that th decreas in O, has a very small 





effect on the default Nusselt number similar to the small 


effect of o,- 
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peo 
Number | X/D | S/D | Changed aa ere rence 
| 4.83 _| 
3980 1 5 25 C2 32.90 32.28 1.92 
3980 1.5 2.5 | Sigma-E S 33.86 32.28 4.89 
3980 1.5 2.5 | Sigma-K -30 34.32 32.28 6.32 
| 7310 | 1.5 | 25 | Default | 0 | 5058 | 4961 | 1.96 | 
7310 1.5 2.5 C2 49.83 49.61 0.44 
7310 1.5 2.5 | Sigma-E = 49.82 49.61 0.42 
7310 1.5 2.5 | Sigma-K 49.95 49.61 0.69 
Freee [is [e8 [Beit [0 rn [rns 
1. 5 2. 5 C2 73, 59 77.75 
1.5 2.5 71.29 77.75 
1.5 2.5 71.17 77.75 


Peso [as [es | Det [0 [esse [e573 fines 
C2 


3590 25 25 23. 89 26.73 


3590 2.5 2.5 29.06 26.73 
3590 2.5 2.5 29.39 26.73 


reo ['25 [as | Deiat [0 [aso [aaa 
C2 


7340 25 25 40. 00 44.99 
7340 2.5 2.5 44.88 44.99 
7340 2.5 2.5 44.88 44.99 


56.06 71 61 

60.89 71.61 

60.44 71.61 
Table 5. Summary of SKE performance when 





and 6, are decreased by 30% 





For the SKE model, general trends observed were that 
for the lower Reynolds numbers (3980, 7310, 3590, 7340) the 
default SKE model overestimated the Nusselt number for the 
array. For higher Reynolds numbers (13800, 13900) the 


trend was reversed and the default SKE model consistently 





underestimated th xperimental Nusselt number. This trend 


seems to be independent of the longitudinal spacing (xX) 





between pins. 


2. The RNG Model 


The next model analyzed was the RNG model. For each 
Reynolds number a run was carried out for the default model 


settings. The resulting Nu was compared to the same 
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experimental Nu from the SKE runs. There are 10 model 
coefficients in the RNG model that can be modified using 


ANSYS 6.0. Only two of these parameters (n, and B) have a 





different ffect on the two transport equations (1.1-2) 





compared to the SKE model. Both of these parameters are 


the multipliers in the RNG formulation for C,, equation 








(1.9). The constant Nn, has the greatest effect of the two 
on the C, formulation. Hence this was selected as the only 


parameter to vary in the RNG trials; the others were 
unchanged from the SKE model since they had been accounted 


for in: those trial runs. <The value of 1). in the RNG model 





was also changed by 30% up and down from its default 


setting. Table 6 summarizes the resulting changes. 


fmt o[ so | Gomes | Ques | "nse | toe 
Number | X/D | S/D Changed eraee Difference 
Default 
Eta-inf 
Eta-inf 
Default 
Eta-inf 
Eta-inf 
Default 
Eta-inf 
Eta-inf 
Default 
Eta-inf 
Eta-inf 
Default 
Eta-inf 
Eta-inf 
Default 
Eta-inf 
Eta-inf 


= 
1 
1 
1 
1 
1 
1 
1 
2 


av th. OF GV ar Oi eh Er OF Oi 





NMNNYUY ND 
ananw#»nuoan o ol 


Table: ©. Summary of RNG Performance When 7, is 


Varied 
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The RNG results for the default model underestimate 


the experimental Nusselt for all but the lowest 2 Reynolds 





numbers. This trend is independent of the X/D ratio. For 


the lowest Reynolds numbers the RNG default model 





overestimates the Nusselt number. The RNG is a more 














sophisticated model in that it takes into account the 
variation of the turbulent kinetic energy in the flow field 
when determining the C, term. However this added parameter 
provides a more accurate Nusselt number than the SKE model 


only at Reynolds numbers less than or equal to 7340. This 


agrees with prior work described by Shah et al (2001) for 














2-D cases using the RNG model in similar flow scenarios. 


Varying the constant NN, does not have predictable results. 





Increasing the constant by 30% always pushes the simulated 








Nusselt number higher. Decreasing the constant has less 








predictable results, but in the extended configuration (X/D 





= 2.5), decreasing Nn, decreased the Nusselt number. This 
appears to be geometry dependent. The results of 
decreasing ‘...-can. be seen im Che change: of the distribution 
of turbulent kinetic energy in the flow field for the 3980 
Reynolds number in Figure 7. Figure 8 shows the difference 
in the temperature field for the same conditions. Both 


views are looking down through the adiabatic plane. 
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MODAL SOLUTION 


Oe=1080,2/0=2.5,%/T=1. 


MODAL SOLUTION 








7H 10 she 


121: 


3-303 
207 413 


5, Tw =306, h/v/p=5/17/2, COLZ, ADS 


7H 10 202 
tL 


2.061 3.428 
1,561 2,601 3,642 


Oe=1080,2/0=2.5,2/TH=1.5, Te =306, h/v/p=5/17/2 , COLe, 250, DNC, etaint-60 


Figure 7. 


Change in the turbulent kinetic energy 


distribution of an RNG run at Re=3980 where nN, 


was decreased by 60% 
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AN. 
7H 16 Moz 
SLs25555 


299.366 300,057 301,769 
299.217 


Re=3980,5/162.5,x/061.5, Tw =306, h/w/pr5/17/2, coms, PRG 


AN. 
7H 16 Moe 
21:41:21 


298.293 300,095 301,506 303.713 
295.196 301,003 


305.52 
304.616 
Re=3980,5/162.5,x/061.5, 1 =306, h/w/pr5/17/2, come, 250, RBs, eteinf-60 





Figure 8. Change in the temperature distribution 
of an RNG run at Re=3980 where N, was decreased 
by 60% 





In other runs it was observed that any change up or 


down in the C, constant brought a change in the Nusselt 


number in the same direction. The constant Cc, was a 





consistent predictor but its effectiveness decreased the 





greater the change in the constant. It was unable to 


change the Nusselt numbers by enough to bring them closer 
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to th xperimental values without causing the transport 


equations to become unstable. 


Other runs also included runs where C, was decreased 





by 30%. The Nusselt number changed in the same direction 


for all Reynolds numbers for a given spacing configuration. 





When X/D=1.5, decreases in C, brought increases in the 


Nusselt number with respect to the default RNG runs. 





However when X/D=2.5, decreases in C, resulted in decreases 


in the Nusselt number. 


3. The NKE Model 


The next model run was the New k-€ (NKE) Model. This 





model is proposed in the literature as being suitable for 





high Reynolds numbers. The model bears this out for its 





default case. The differenc between th simulated and 


experimental Nusselt numbers decreases as Reynolds number 





increases but only for the cases of X/D=1.5. At Re=13900 


the numerical runs were unsuccessful. The results are 
summarized in Table 7. The values of Cy were chosen to be 
evenly spaced between 0.43 and 1.0. The lower number, 


0.43, corresponds with the default setting, and 1.0 


corresponds with the limit of = as TN] goes to infinity. 
H+ 
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Reynolds 
Number | X/ 


3980 
3980 
3980 
7310 
7310 
7310 


Variable Simulated | Experimental | Percent 
S/D | Changed Value Nu Nu Diff. 


Default 0.43 35.94 32.28 
CiM 0.6 31.44 32.28 
CiM 0.8 32.88 32.28 

Default 0.43 53.53 49.61 
CiM 0.6 46.91 49.61 
CiM 0.8 50.71 49.61 

Default 0.43 74.87 77.75 
CiM 0.6 70.23 77.75 
CiM 0.8 75.97 77.75 

Default 0.43 31.66 26.73 
C1iM 0.6 25.17 26.73 
C1iM 0.8 25.00 26.73 

Default 0.43 46.45 44.99 
CiM 0.6 38.67 44.99 
CiM 0.8 39.07 44.99 

Default 0.43 71.61 
CiM 0.6 71.61 
CiM 0.8 71.61 


iw) 
Ny ~ Py wy Y|y yw 
aj a ayo oO ayo oOo on 


1.5 
1.5 
1.5 
1.5 
1.5 
1.5 
1.5 
1.5 
1.5 





Table 7. Summary of NKE Performance 





Changing the C,, constant for the X/D ‘case of 1.5 


results in improved performance between 0.6 and 0.8, 





increasing to 0.8 for the highest Reynolds number of the 











K/D. = “Les, Case. In. the. AAD. = 20. Case: changing: the -€.. 


constant has a mixed effect. 
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Vv. CONCLUSIONS AND RECOMENDATIONS 


A. OPTIMIZATION OF THE SHORT PIN FIN HEAT EXCHANGER 


Table 8 shows the optimal constants for the short pin 





fin heat exchanger. There is substantial improvement in 
the X/D=1.5 spacing ratio model. The lower Reynolds number 
flows, regardless of spacing ratio, are not well 


represented by the SKE default model. The SKE default model 


consistently overestimates the amount of turbulent kinetic 











energy in low Reynolds number flows resulting in simulated 


flows with high Nusselt numbers. The ineréase in: the C;. 


constant and the reduction in the C, constant weight terms 


in their specific models to decrease turbulent kinetic 


energy and increase the dissipation rate. This can also be 





seen in Table 4 where C is increased by 30% and this 


le 
drives the dissipation rate up as the transport equations 


(1.1-.2) are solved resulting in a lower Nusselt number, 








closer to the experimental value. 


Variable | Simulated |Experimental} SKE Percent 
Reynolds Combination Nu Nu Def. Nu_ {Improvement 


| 7340 | SKE | Defautt_ | 45.69 | 4499 | 45.69 | 0.00 





Table 8. Optimal Constants for the Cross Flow 
Pin Fin Heat Exchanger. 
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Table 9 shows the performance of the optimal constants 
against Metzger's values. 
Reynolds Variable | Simulated |Experimental} Percent 
Number Combination Nu Nu Difference 
3980 C1430 32.19 32.28 
7310 Default 49.40 49.61 


SKE 
SKE 
13900_| RNG 





Table 9. Optimal Constants Compared with 
Metzger's Experimental Values 





B. FUTURE PLANS 


The next steps for this study should be to explore 
changes to terms in the models that affect the balance of 


turbulent kinetic energy and the turbulent kinetic energy 








dissipation rate (1.1-.2). This balance seems to have the 
most significant effect on the performance of a turbulence 
model. Specific changes need to be in exploring a reduction 


of the C, constant for higher Reynolds numbers in the 





expanded configuration (X/D=2.5) of the SKE model. There is 








a possibility that the viscous dissipation term becomes 
greater as the flow moves farther downstream and away from 
turbulence creating obstructions such as the pins. Table 4 


shows that the changes of +/- 30% of C bracket the 


le 





experimental Nusselt number with the exception of the 





Re=7310 run. This leads to varying constants in the model 





that are directly tied to source and dissipation terms as 











Ci. 282 Iterating the percent change should bring that 





constant into close agreement with the experimental value. 


Another specific area is the higher Reynolds number 
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performance of the RNG model. The N, constant does not 


provide predictable improvements for the model. At higher 





Reynolds numbers there may be a benefit in increasing the 


C, term that reduces the contribution of the viscous 
dissipation rate, which appears as ¢’. 


For future considerations, temperature differences 


must be increased to more accurately model environmental 





temperatures, H/D ratios should be varied slightly (0.5- 





2.0) to s where these relations need to be improved. 


Another environmental concern is to model the 





rotational environment experienced in a gas turbine's power 





turbine. For further heat transfer performance refinement, 
attention will have to be paid to the variations available 


for the wall turbulence model. 


Many compact heat exchangers use different pin shapes 
than cylinders. The ANSYS code has been written with 
vertical planes midway between each row that will 
facilitate changing entire rows to different span wise 


spacing and pin shapes. 
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